MODULE constants
  IMPLICIT NONE
  REAL*8,PARAMETER:: c_pi=3.14159265358979323846d0
  !REAL*8,PARAMETER:: c_pi=4.d0*datan(1.d0)
  REAL*8,PARAMETER:: c_try=1.23456789123456789123456d0
  REAL*8,PARAMETER:: c_na = 6.02214d+23 ! Avogadro's number
  REAL*8,PARAMETER:: c_h = 6.62606896d-34 ! Planck's constant (J*s)
  REAL*8,PARAMETER:: c_c = 2.99792458d+10 ! Speed of light in vacuum (cm/s)
  REAL*8,PARAMETER:: c_h2wn = 2.194746d+05 ! Hartree to cm^-1
  REAL*8,PARAMETER:: c_wn2h = 1.d0/2.194746d+05 ! cm^-1 to hartree
  REAL*8,PARAMETER:: c_amu = 1.660538782d-27 ! = 1/_na/1000
  REAL*8,PARAMETER:: c_gamma = 4.d0 * c_pi * c_pi * c_c * 1.0d-23 / c_h / c_na;
  REAL*8,PARAMETER:: c_bohr=5.291772108d-11!bohr radius in meters =
  REAL*8,PARAMETER:: c_angs2bohr=1d-10/c_bohr
  REAL*8,PARAMETER:: c_emass=9.10938215d-31!electron mass in kg
  REAL*8,PARAMETER:: c_amu2emass=c_amu/c_emass
END MODULE constants

MODULE global
  IMPLICIT NONE
  INTEGER maxcell,maxk,Ndof,Ndim,Natom,Ncell,N_k,nMR,vciorder
  INTEGER, ALLOCATABLE::maxbasis(:)
  REAL*8, ALLOCATABLE ::xmass(:),k_freq(:,:),hrmfreq(:),NC(:,:),freq(:)
  COMPLEX*16, ALLOCATABLE ::k_coef(:,:,:),cck_coef(:,:,:)
  REAL*8  :: geo(100)
  LOGICAL::vscf,vci,vpt,vav,summary
END MODULE global

MODULE forceconstants
  IMPLICIT NONE

  REAL*8,ALLOCATABLE                              :: x_grad(:)
  COMPLEX*16,ALLOCATABLE                        :: q_grad(:,:)
  REAL*8,ALLOCATABLE                          :: x_harm(:,:,:)
  COMPLEX*16,ALLOCATABLE                    :: q_harm(:,:,:,:)
  !  real*8,allocatable                      :: x_cube(:,:,:,:,:)
  !  complex*16,allocatable               :: q_cube(:,:,:,:,:,:)
  !  real*8,allocatable                 :: x_quar(:,:,:,:,:,:,:)
  !  complex*16,allocatable :: q_quar(:,:,:,:,:,:,:,:)

  !variables for nMR QFF, n=1,2,3,4
  CHARACTER*100 :: title1MR,title2MR,title3MR
  REAL*8:: E0
  REAL*8,ALLOCATABLE :: Gi(:),Hii(:),Ciii(:),Qiiii(:)
  REAL*8,ALLOCATABLE :: Hij(:,:),Ciij(:,:),Qiijj(:,:),Qiiij(:,:)
  REAL*8,ALLOCATABLE :: Cijk(:,:,:),Qiijk(:,:,:)
  REAL*8,ALLOCATABLE :: Qijkm(:,:,:,:)
  ! real*8 :: Gi(Ndof),Hii(Ndof),Ciii(Ndof),Qiiii(Ndof)
  ! real*8 :: Hij(Ndof,Ndof),Ciij(Ndof,Ndof),Qiijj(Ndof,Ndof),Qiiij(Ndof,Ndof)
  ! real*8 :: Cijk(Ndof,Ndof,Ndof),Qiijk(Ndof,Ndof,Ndof)
  ! real*8 :: Qijkm(Ndof,Ndof,Ndof,Ndof)

END MODULE forceconstants

MODULE integral
  IMPLICIT NONE
CONTAINS

  REAL*8 FUNCTION q_integral(level,mode,power,diff)
    ! <n|Q_m^power|n-diff>
    USE constants
    USE Global
    IMPLICIT NONE
    INTEGER ::level,mode,power,diff
    REAL*8::n,n_2,n_3
    REAL*8::gamma,gamma_2,gamma_3
    ! <n|Q_m^power|n-diff>
    n=level
    n_2 = n * n;
    n_3 = n * n * n;
    gamma = freq(mode)
    gamma_2 = gamma * gamma
    gamma_3 = gamma * gamma * gamma
    q_integral=0.d0

    SELECT CASE (power)

    CASE (1)
       IF (diff == 1) THEN
          q_integral=DSQRT(n * 0.5d0 / gamma);
       END IF
    CASE (2)
       IF (diff == 0) THEN
          q_integral=(n + 0.5d0) / gamma;
       ELSE IF (diff == 2)THEN
          q_integral=0.5d0 * DSQRT(n * (n - 1.d0)) / gamma;
       ENDIF
    CASE (3)
       IF (diff == 1)THEN
          q_integral=3.d0 * DSQRT(n_3 / (8.d0 * gamma_3));
       ELSE IF (diff == 3) THEN
          q_integral=DSQRT(0.125d0 * n * (n - 1.d0) * (n - 2.d0) / gamma_3);
       ENDIF
    CASE (4)
       IF (diff == 0)THEN
          q_integral=0.75d0 * (n_2 * 2.d0 + n * 2.d0 + 1.d0) / gamma_2;
       ELSE IF (diff == 2)THEN
          q_integral=0.5d0 * (2.d0 * n - 1.d0) * DSQRT(n * (n - 1.d0)) / gamma_2;
       ELSE IF (diff == 4)THEN
          q_integral=0.25d0 * dsqrt(n * (n - 1.d0) * (n - 2.d0) * (n - 3.d0)) / gamma_2;
       ENDIF
    END SELECT
  END FUNCTION q_integral

  REAL*8 FUNCTION kin_integral(level,m,diff)
    USE constants,ONLY:c_gamma
    USE Global
    IMPLICIT NONE
    INTEGER ::m,diff,level,n,n_2,n_3
    REAL*8::gamma,gamma_2,gamma_3
    ! <n|del^2/del(q)^2|n-diff>
    n=level
    n_2 = n * n;
    n_3 = n * n * n;
    gamma = c_gamma * hrmfreq(m)
    gamma_2 = gamma * gamma
    gamma_3 = gamma * gamma * gamma
    IF (diff == 0)THEN
       kin_integral=0.5 * gamma * (n + 0.5)
    ELSE IF (diff == 2)THEN
       kin_integral=0.25 * gamma * SQRT((n + 1.0)*dfloat(n) )
    ELSE
       kin_integral=0.0
    ENDIF
  END FUNCTION kin_integral
END MODULE integral

MODULE integral2
  USE integral
  IMPLICIT NONE
CONTAINS

  REAL*8 FUNCTION fullIntegral(mode,power,coef)
    USE global
    IMPLICIT NONE
    INTEGER ::mode,power,n1,n2
    REAL*8::coef(nDOF,maxbasis(mode))
    fullIntegral = 0.d0;
    DO n1 = 0,maxbasis(mode)-1
       DO n2 = 0,maxbasis(mode)-1
          fullintegral = fullintegral + coef(mode,n1+1) * coef(mode,n2+1) * q_integral(MAX(n1, n2), mode, power, ABS(n2-n1))
       ENDDO
    ENDDO
  END FUNCTION fullintegral

END MODULE integral2

MODULE qVSCFintegrals
  USE global,ONLY:ndof,freq
  USE forceconstants
  USE integral
  IMPLICIT NONE
CONTAINS
  REAL*8 FUNCTION uzero(state,mode)
    IMPLICIT NONE
    INTEGER::m1,m2,mode,state(nDOF)
    uzero=0.d0
    DO m1=1,ndof
       IF (m1 .EQ. mode)CYCLE
       !print*,'m1= ',m1
       !uzero= uzero + 0.5*Fii(m1) * q_integral(state(m1), m1, 2, 0)
       uzero= uzero + hii(m1) * q_integral(state(m1), m1, 2, 0)
       !print*,'uzero= ',uzero
       DO m2=1,ndof
          IF (m2 .EQ. mode .OR. m1 .EQ. m2)CYCLE
          !uzero= uzero + 0.125d0 * Fiijj(m1,m2) * q_integral(state(m1), m1, 2, 0)* q_integral(state(m2), m2, 2, 0)
          uzero= uzero + 0.5d0 * qiijj(m1,m2) * q_integral(state(m1), m1, 2, 0)* q_integral(state(m2), m2, 2, 0)
       ENDDO
    ENDDO
  END FUNCTION uzero
  REAL*8 FUNCTION utwo(state,mode)
    IMPLICIT NONE
    INTEGER::m,mode,state(nDOF)
    utwo=0.d0
    DO m=1,ndof
       IF (m .EQ. mode)CYCLE
       !utwo= utwo + 0.5d0 *Fiijj(mode,m) * q_integral(state(m), m, 2, 0)
       utwo= utwo + 2.0d0 *qiijj(mode,m) * q_integral(state(m), m, 2, 0)
    ENDDO
    !utwo= utwo + Fii(mode)
    utwo= utwo + 2.d0*hii(mode)
  END FUNCTION utwo

  REAL*8 FUNCTION Ukinetic(state,mode)
    IMPLICIT NONE
    INTEGER ::m,state(nDOF),mode
    ! <n|del^2/del(q)^2|n-diff>
    Ukinetic=0.d0
    DO m=1,ndof
       IF (m .EQ. mode)CYCLE
       Ukinetic= Ukinetic+0.5 * freq(m)* (state(m) + 0.5)
    ENDDO
    RETURN
  END FUNCTION Ukinetic

  REAL*8 FUNCTION Ekinetic(state)
    IMPLICIT NONE
    INTEGER ::m,state(nDOF)
    ! <n|del^2/del(q)^2|n-diff>
    Ekinetic=0.d0
    DO m=1,ndof
       Ekinetic= Ekinetic+0.5 * freq(m)* (state(m) + 0.5)
    ENDDO
    RETURN
  END FUNCTION Ekinetic

  REAL*8 FUNCTION Epotential(state)
    IMPLICIT NONE
    INTEGER::m1,m2,state(nDOF)
    Epotential=0.d0
    DO m1=1,ndof
       !Epotential= Epotential + 0.5*Fii(m1) * q_integral(state(m1), m1, 2, 0)
       Epotential= Epotential + hii(m1) * q_integral(state(m1), m1, 2, 0)
       DO m2=1,ndof
          IF ( m1 .EQ. m2)CYCLE
          !Epotential= Epotential + 0.125d0 * Fiijj(m1,m2) * q_integral(state(m1), m1, 2, 0)* q_integral(state(m2), m2, 2, 0)
          Epotential= Epotential + 0.5d0 * qiijj(m1,m2) * q_integral(state(m1), m1, 2, 0)* q_integral(state(m2), m2, 2, 0)
       ENDDO
    ENDDO
  END FUNCTION Epotential

END MODULE qVSCFintegrals

MODULE yazar
  IMPLICIT NONE
CONTAINS
  SUBROUTINE writematrix(order,matrix)
    IMPLICIT NONE
    INTEGER i,j,order
    REAL*8,DIMENSION(order,order)::matrix
    DO i=1,order
       WRITE(*,'(10D20.9,6x)',ADVANCE='no')(matrix(i,j),j=1,order)
       WRITE(*,*)
    ENDDO
  END SUBROUTINE writematrix
END MODULE yazar


MODULE modvci
CONTAINS
  SUBROUTINE statelabel(nthstate,ketlabel)
    USE global
    IMPLICIT NONE
    INTEGER::m,state,nthstate,ketlabel(ndof)
    state = nthState;
    DO m=1,ndof
       ketlabel(m) = MOD(state,maxBasis(m));
       state = state / maxBasis(m);
    ENDDO
    RETURN
  END SUBROUTINE statelabel

  SUBROUTINE fundiffer(bralabel,ketlabel,differvec,ndiffer)
    USE global
    IMPLICIT NONE
    INTEGER::i,m, bralabel(ndof),ketlabel(ndof),differvec(ndof),ndiffer
    i = 0;
    DO m=1,ndof
       IF (braLabel(m) /= ketLabel(m))THEN
          differvec(i+1) = m;
          i=i+1;
       ENDIF
    ENDDO
    ndiffer=i
  END SUBROUTINE fundiffer
END MODULE modvci

MODULE lapack
  IMPLICIT NONE
CONTAINS

  SUBROUTINE diag(order,matrix,eigenvalues)
    IMPLICIT NONE
    INTEGER :: order,info,lda
    REAL(8), DIMENSION(order,order)::matrix
    REAL(8), DIMENSION(order)::eigenvalues
    REAL(8), DIMENSION(3*order-1) :: work
    lda=3*order-1
    CALL dsyev('V','U',order,matrix,order,eigenvalues,work,lda,info)
    !              PRINT*,'info=',info
  END SUBROUTINE diag

  SUBROUTINE complex_diag(order,matrix,eigenvalues)
    IMPLICIT NONE
    INTEGER :: order,info,lda,Lwork
    COMPLEX*16, DIMENSION(order,order)::matrix
    REAL(8), DIMENSION(order)::eigenvalues
    COMPLEX*16, DIMENSION(2*order-1) :: work
    REAL(8), DIMENSION(3*order-2) :: Rwork
    lda=3*order-1
    Lwork=2*order-1
    CALL zheev('V','U',order,matrix,order,eigenvalues,work,Lwork,Rwork,info)
  END SUBROUTINE complex_diag

END MODULE lapack

!*********subroutines**************************************************************************************

SUBROUTINE myvscf(state,energy,iter)
  USE global
  USE lapack
  USE yazar
  USE constants
  IMPLICIT NONE
  REAL*8:: ham(MAXVAL(maxbasis),MAXVAL(maxbasis)),energies(MAXVAL(maxbasis)),coef(nDOF,MAXVAL(maxbasis))
  INTEGER:: mode,iter,maxiter,state(nDOF)
  REAL*8:: energy,oldenergy
  REAL*8,PARAMETER:: thresh=0.000001
  iter=0
  maxiter=20
  oldenergy=1000
  energy=0.0
  coef=0.0
  !Forming initial wavefunction using HO basis.
  DO mode=1,nDOF
     coef(mode,state(mode)+1)=1.d0
  ENDDO
  DO WHILE ((ABS(energy-oldenergy) .GT. thresh) .AND. iter .LT. maxiter )
     oldenergy=energy
     energy=0.0
     !   PRINT*,'**************************  iteration **************************',iter
     iter=iter+1
     DO mode=1,nDOF
        CALL meanFieldHamiltonian(mode,coef,ham)
        !   PRINT*,'hamiltonian for mode ',mode
      !  CALL writematrix(maxbasis(mode),ham*c_h2wn)
        CALL diag(maxbasis(mode),ham,energies)
        energy=energies(state(mode)+1)
        ! PRINT*,"energy= ", energies*c_h2wn!,energies(1)
        coef(mode,:)=ham(:,state(mode)+1)
        ! PRINT*,'eigvec'!ham(1,:)
     ENDDO
  ENDDO
  ! print*,'vscf done for', state
  RETURN
END SUBROUTINE myvscf

SUBROUTINE qvscf(state,energy,iter)
  USE global
  USE lapack
  USE yazar
  USE constants
  USE qVSCFintegrals
  IMPLICIT NONE
  INTEGER:: mode,iter,maxiter,state(nDOF)
  REAL*8:: energy,oldenergy
  REAL*8,PARAMETER:: thresh=0.000001
  iter=0
  maxiter=20
  oldenergy=1000
  energy=0.0
  !freq=hrmfreq
  DO WHILE ((ABS(energy-oldenergy) .GT. thresh) .AND. iter .LT. maxiter )
     oldenergy=energy
     iter=iter+1
     DO mode=1,nDOF
        freq(mode)=dsqrt(utwo(state,mode))
        ! print*,'freq',freq(mode)
     ENDDO
     energy=Ekinetic(state)+Epotential(state)
  ENDDO

END SUBROUTINE qvscf

SUBROUTINE meanFieldHamiltonian(mode,coef,Hvscf)
  USE global
  USE forceconstants
  USE integral
  USE integral2
  IMPLICIT NONE
  INTEGER::mode,diff,i,j,m,m1,m2
  REAL*8::tmp,meanfieldpotential,meanpotkin
  REAL*8::Hvscf(maxbasis(mode),maxbasis(mode)),coef(nDOF,maxbasis(mode))
  !PRINT*, "forming VSCF Hamiltonian for mode" , mode
  meanpotkin=meanFieldPotential(mode,coef)
  ! PRINT*,'meanpotkin',meanpotkin
  Hvscf=0.d0
  DO i=0, maxbasis(mode)-1
     tmp=0.0
     tmp= tmp + (i + 0.5) * hrmfreq(mode) * 0.5;!kinetic term
     tmp= tmp + Hii(mode) * q_integral(i, mode, 2, 0);
     tmp= tmp + Qiiii(mode) * q_integral(i, mode, 4, 0);
     tmp= tmp + meanpotkin


     IF (nMR .GT. 1) THEN
        DO m=1,ndof
           IF (m .EQ. mode)CYCLE
           tmp= tmp + hij(m,mode) * q_integral(i, mode, 1, 0) * fullintegral(m, 1,coef);
           tmp= tmp + qiijj(m,mode) * q_integral(i, mode, 2, 0) * fullintegral(m, 2,coef);
           tmp= tmp + ciij(mode,m) * q_integral(i, mode, 2, 0) * fullintegral(m, 1,coef);
           tmp= tmp + ciij(m,mode) * q_integral(i, mode, 1, 0) * fullintegral(m, 2,coef);
           tmp= tmp + qiiij(mode,m) * q_integral(i, mode, 3, 0) * fullintegral(m, 1,coef);
           tmp= tmp + qiiij(m,mode) * q_integral(i, mode, 1, 0) * fullintegral(m, 3,coef);
        ENDDO
        IF (nMR .GT. 2) THEN
           DO m1=1,ndof-2
              IF (m1 .EQ. mode)CYCLE
              DO m2=m1+1,ndof-1
                 IF (m2 .EQ. mode) CYCLE
                 tmp= tmp + cijk(mode,m1,m2) * q_integral(i, mode, 1, 0) * fullintegral(m1, 1,coef)  * fullintegral(m2, 1,coef);
                 tmp= tmp + qiijk(mode,m1,m2) * q_integral(i, mode, 2, 0) * fullintegral(m1, 1,coef) * fullintegral(m2, 1,coef);
                 tmp= tmp + qiijk(m1,mode,m2) * q_integral(i, mode, 1, 0) * fullintegral(m1, 2,coef)* fullintegral(m2, 1,coef);
                 tmp= tmp + qiijk(m2,mode,m1) * q_integral(i, mode, 1, 0) * fullintegral(m1, 1,coef) * fullintegral(m2, 2,coef);
              ENDDO!!m2
           ENDDO!m1
        ENDIF!nMR>2
     ENDIF!nMR>1
     Hvscf(i+1,i+1)=tmp
     DO j=i+1,maxbasis(mode)-1!
        tmp=0.0
        diff = j - i;
        IF (diff > 4) EXIT !break;
        IF (diff .EQ. 2) tmp= tmp-SQRT(float(j * (j - 1))) * hrmfreq(mode) * 0.25 !!kinetic term
        tmp= tmp + q_integral(j, mode, 1, diff) * gi(mode);
        tmp= tmp + q_integral(j, mode, 2, diff) * hii(mode);
        tmp= tmp + q_integral(j, mode, 3, diff) * ciii(mode);
        tmp= tmp + q_integral(j, mode, 4, diff) * qiiii(mode);
        ! print*,'test2',tmp
        IF (nMR .GT. 1)THEN
           DO m=1,ndof
              IF (m .EQ. mode)CYCLE
              tmp= tmp + hij(m,mode) * q_integral(j, mode, 1, diff) * fullintegral(m, 1,coef);
              tmp= tmp + qiijj(m,mode) * q_integral(j, mode, 2, diff) * fullintegral(m, 2,coef);
              tmp= tmp + ciij(mode,m) * q_integral(j, mode, 2, diff) * fullintegral(m, 1,coef);
              tmp= tmp + ciij(m,mode) * q_integral(j, mode, 1, diff) * fullintegral(m, 2,coef);
              tmp= tmp + qiiij(mode,m) * q_integral(j, mode, 3, diff) * fullintegral(m, 1,coef);
              tmp= tmp + qiiij(m,mode) * q_integral(j, mode, 1, diff) * fullintegral(m, 3,coef);
           ENDDO
           IF (nMR .GT. 2) THEN
              DO m1=1,ndof
                 IF (m1 .EQ. mode)CYCLE
                 DO m2=m1 + 1,ndof
                    IF (m2 .EQ. mode)CYCLE
                    tmp= tmp + cijk(mode,m1,m2)*q_integral(j,mode,1,diff) * fullintegral(m1,1,coef)* fullintegral(m2,1,coef);
                    tmp= tmp + qiijk(mode,m1,m2)*q_integral(j,mode,2,diff) * fullintegral(m1,1,coef)* fullintegral(m2,1,coef);
                    tmp= tmp + qiijk(m1,mode,m2)*q_integral(j,mode,1,diff) * fullintegral(m1,2,coef)* fullintegral(m2,1,coef);
                    tmp= tmp + qiijk(m2,mode,m1)*q_integral(j,mode,1,diff) * fullintegral(m1,1,coef)* fullintegral(m2,2,coef);
                 ENDDO!!m2
              ENDDO!!m1
           ENDIF!nMR>2
        ENDIF!nMR>1
        ! print*,'test3',tmp
        Hvscf(i+1,j+1)=tmp
     ENDDO!!end for j
  ENDDO!!end for i
END SUBROUTINE meanFieldHamiltonian

REAL*8 FUNCTION meanfieldpotential(mode,coef)
  USE global
  USE forceconstants
  USE integral
  USE integral2
  USE constants
  IMPLICIT NONE
  REAL*8 :: tmp,meanFieldKinetic,coef(nDOF,maxbasis(mode))
  INTEGER::m,m1,m2,m3,mode
  tmp = 0.0;
  ! WRITE(*,'(d50.25)')c_gamma
  DO m=1,Ndof
     IF (m .EQ. mode)CYCLE
     tmp =tmp+ fullintegral(m, 1,coef) * gi(m);
     tmp =tmp+ fullintegral(m, 2,coef) * hii(m);
     tmp =tmp+ fullintegral(m, 3,coef) * ciii(m);
     tmp =tmp+ fullintegral(m, 4,coef) * qiiii(m);
     tmp =tmp+ meanFieldKinetic(m,coef);
  ENDDO
  IF (nMR .GT. 1)THEN
     ! print*, "**************************nMR >1 "
     DO m1=1,nDOF
        IF (m1 .EQ. mode)CYCLE
        DO m2=m1+1,nDOF
           IF (m2 .EQ. mode .OR. m1 .EQ. m2)CYCLE
           tmp =tmp+ fullintegral(m1, 1,coef) * fullintegral(m2, 1,coef) * hij(m1,m2);
           tmp =tmp+ fullintegral(m1, 2,coef) * fullintegral(m2, 2,coef) * qiijj(m1,m2);
           tmp =tmp+ fullintegral(m1, 2,coef) * fullintegral(m2, 1,coef) * ciij(m1,m2);
           tmp =tmp+ fullintegral(m1, 1,coef) * fullintegral(m2, 2,coef) * ciij(m2,m1);
           tmp =tmp+ fullintegral(m1, 3,coef) * fullintegral(m2, 1,coef) * qiiij(m1,m2);
           tmp =tmp+ fullintegral(m1, 1,coef) * fullintegral(m2, 3,coef) * qiiij(m2,m1);
        ENDDO!m2
     ENDDO!m1
     IF (nMR .GT. 2)THEN
        DO m1=1,nDOF
           IF (m1 .EQ. mode) CYCLE
           DO m2=m1+1,nDOF
              IF (m2 .EQ. mode) CYCLE
              DO m3=m2+1,nDOF
                 IF (m3 .EQ. mode) CYCLE
                 tmp =tmp+ fullintegral(m1, 1,coef) * fullintegral(m2, 1,coef) * fullintegral(m3, 1,coef) * cijk(m1,m2,m3);
                 tmp =tmp+ fullintegral(m1, 2,coef) * fullintegral(m2, 1,coef) * fullintegral(m3, 1,coef) * qiijk(m1,m2,m3);
                 tmp =tmp+ fullintegral(m1, 1,coef) * fullintegral(m2, 2,coef) * fullintegral(m3, 1,coef) * qiijk(m2,m1,m3);
                 tmp =tmp+ fullintegral(m1, 1,coef) * fullintegral(m2, 1,coef) * fullintegral(m3, 2,coef) * qiijk(m3,m1,m2);
              ENDDO!m3
           ENDDO!m2
        ENDDO!m1
     ENDIF!nMR>2
  ENDIF!nMR>1
  meanfieldpotential=tmp
  !  print*, tmp*c_h2wn,"meanfield pot for" ,mode
  RETURN
END FUNCTION meanfieldpotential

REAL*8 FUNCTION meanfieldkinetic(mode,coef)
  USE global
  USE forceconstants
  USE integral
  USE integral2
  USE constants
  IMPLICIT NONE
  REAL*8 :: tmp,coef(nDOF,maxbasis(mode))
  INTEGER::i,j,diff,mode

  tmp = 0.0;
  DO i=0, maxbasis(mode)-1
     DO j=0, maxbasis(mode)-1
        diff = j - i;
        IF (diff .EQ. 0) THEN
           tmp =tmp+ coef(mode,i+1) * coef(mode,i+1) * (i + 0.5) * hrmfreq(mode) * 0.5;
        ELSE IF (ABS(diff) .EQ. 2) THEN
           tmp =tmp+ coef(mode,i+1) * coef(mode,j+1) * (-SQRT(float((MAX(i, j)) * (MAX(i, j) - 1)))) * hrmfreq(mode) * 0.25;
        ENDIF
     ENDDO!j
  ENDDO!i
  meanfieldkinetic=tmp
  ! print*, tmp*c_h2wn,"meanfield kinetic for" ,mode
  RETURN
END FUNCTION meanfieldkinetic

SUBROUTINE init_QFF
  USE global, ONLY : nDOF,nMR
  USE forceconstants
  IMPLICIT NONE
  ALLOCATE(Gi(Ndof),Hii(ndof),Ciii(ndof),Qiiii(ndof))
  Gi=0.d0;Hii=0.d0;Ciii=0.d0;Qiiii=0.d0
  IF (nMR .GT. 1) THEN
     ALLOCATE(Hij(ndof,ndof),Ciij(ndof,ndof),Qiijj(ndof,ndof),Qiiij(ndof,ndof))
     Hij=0;Ciij=0;Qiijj=0;Qiiij=0
     IF (nMR .GT. 2) THEN
        ALLOCATE(Cijk(ndof,ndof,ndof),Qiijk(ndof,ndof,ndof))
        Cijk=0;Qiijk=0
     ENDIF
  ENDIF
END SUBROUTINE init_QFF

SUBROUTINE read_QFF
  USE global, ONLY : geo,nDOF,nMR
  USE forceconstants
  IMPLICIT NONE
  INTEGER :: i,j,k,n,tmp


  OPEN(UNIT=7,FILE='001.hs',STATUS='OLD')
  READ(7,*,END=666)!# Energy / hartree
  READ(7,'(3x,D18.12)',END=666) E0
  READ(7,*,END=666)!# Geometry / Angs amu1/2
  j=MOD(Ndof,3)
  IF(j==0) THEN
     n=Ndof/3
     DO i=1,n
        READ(7,*,END=666) (geo(3*(i-1)+j),j=1,3)
     END DO
  ELSE
     n=(Ndof-j)/3
     DO i=1,n
        READ(7,*,END=666) (geo(3*(i-1)+j),j=1,3)
     END DO
     READ(7,*,END=666) (geo(j),j=3*n+1,Ndof)
  ENDIF

  !   >> >> 1 MR << <<
  READ(7,'(5x,a)') title1MR
  READ(7,*,END=666)!# Gradient / hartree Angs^-1 amu^-1/2
  DO n=1,Ndof
     READ(7,100)i, Gi(n)
  END DO
  READ(7,*,END=666)!# Hessian(i,i) / hartree Angs^-2 amu^-1
  DO n=1,Ndof
     READ(7,100,END=666) i,Hii(n)
  END DO
  READ(7,*,END=666)!# Cubic(i,i,i) / hartree Angs^-3 amu^-3/2
  DO n=1,Ndof
     READ(7,100,END=666)i, Ciii(n)
  END DO
  READ(7,*,END=666)!# Quartic(i,i,i,i) / hartree Angs^-4 amu^-2
  DO n=1,Ndof
     READ(7,100,END=666)i, Qiiii(n)
  END DO
100 FORMAT(i4,D20.10)

  !   >> >> 2 MR << <<
  IF(nMR .GT. 1) THEN
     tmp = ndof * (ndof - 1) / 2

     READ(7,'(5x,a)',END=222) title2MR

     READ(7,*,END=222)!# Hessian(i,j) / hartree Angs^-2 amu^-1
     DO n=1,tmp
        READ(7,200,END=666)i,j, Hij(i,j)
        Hij(j,i)=Hij(i,j)
     END DO

     READ(7,*,END=666)!# Quartic(i,i,j,j) / hartree Angs^-4 amu^-2
     DO n=1,tmp
        READ(7,200,END=666)i,j, Qiijj(i,j)
        Qiijj(j,i)=Qiijj(i,j)
     END DO

     READ(7,*,END=666)!# Cubic(i,i,j) / hartree Angs^-3 amu^-3/2
     DO n=1,2*tmp
        READ(7,200,END=666)i,j, Ciij(i,j)
     END DO

     READ(7,*,END=666)!# Quartic(i,i,i,j) / hartree Angs^-4 amu^-2
     DO n=1,2*tmp
        READ(7,200,END=666)i,j, Qiiij(i,j)
     END DO
200  FORMAT(2i4,D20.10)

     !   >> >> 3 MR << <<
     IF(nMR .GT. 2) THEN
        tmp = ndof * (ndof - 1) * (ndof - 2) / 6

300     FORMAT(3i4,D20.10)
        READ(7,'(5x,a)',END=222) title3MR
        READ(7,*,END=333)!# Cubic(i,j,k) / hartree Angs^-3 amu^-3/2
        DO n=1,tmp
           READ(7,300,END=666)i,j,k, Cijk(i,j,k)
           Cijk(i,k,j)=Cijk(i,j,k)
           Cijk(j,i,k)=Cijk(i,j,k)
           Cijk(j,k,i)=Cijk(i,j,k)
           Cijk(k,i,j)=Cijk(i,j,k)
           Cijk(k,j,i)=Cijk(i,j,k)
        END DO

        READ(7,*)!# Quartic(i,i,j,k) / hartree Angs^-4 amu^-2
        DO n=1,3*tmp
           READ(7,300,END=666)i,j,k, Qiijk(i,j,k)
           Qiijk(i,k,j)=Qiijk(i,j,k)
        END DO
        !  nMR=3
     ENDIF
  ENDIF
  RETURN

222 CONTINUE
  !  nMR=1
  RETURN
333 CONTINUE
  !  nMR=2
  RETURN
666 STOP 'error in QFF.pot'
END SUBROUTINE read_QFF

SUBROUTINE write_QFF
  USE global, ONLY : geo,nDOF
  USE forceconstants
  IMPLICIT NONE
  INTEGER :: i,j,k,n,tmp
  OPEN(UNIT=8,FILE='fc.poly')
  WRITE(8,1000)
  WRITE(8,'(3x,D18.12)') E0
  WRITE(8,1010)!# Geometry / Angs amu1/2
  j=MOD(Ndof,3)
  IF(j==0) THEN
     n=Ndof/3
     DO i=1,n
        WRITE(8,999) (geo(3*(i-1)+j),j=1,3)
     END DO
  ELSE
     n=(Ndof-j)/3
     DO i=1,n
        WRITE(8,999) (geo(3*(i-1)+j),j=1,3)
     END DO
     WRITE(8,999) (geo(j),j=3*n+1,Ndof)
  ENDIF
999 FORMAT(3d18.7)
1000 FORMAT('# Energy / hartree ')
1001 FORMAT('# Dipole / debye ')
1010 FORMAT('# Geometry / Angs amu1/2')

  !   >> >> 1 MR << <<

  WRITE(8,101) title1MR

  WRITE(8,110)!# Gradient / hartree Angs^-1 amu^-1/2
  DO n=1,Ndof
     WRITE(8,100)n, Gi(n)
  END DO

  WRITE(8,120)!# Hessian(i,i) / hartree Angs^-2 amu^-1
  DO n=1,Ndof
     WRITE(8,100)n, Hii(n)
  END DO

  WRITE(8,130)!# Cubic(i,i,i) / hartree Angs^-3 amu^-3/2
  DO n=1,Ndof
     WRITE(8,100)n, Ciii(n)
  END DO

  WRITE(8,140)!# Quartic(i,i,i,i) / hartree Angs^-4 amu^-2
  DO n=1,Ndof
     WRITE(8,100)n, Qiiii(n)
  END DO

100 FORMAT(i4,D20.10)
101 FORMAT('# 1MR ',a)
110 FORMAT('# Gradient / hartree Angs^-1 amu^-1/2 ')
111 FORMAT('# Gradient / debye Angs^-1 amu^-1/2 ')
120 FORMAT('# Hessian(i,i) / hartree Angs^-2 amu^-1 ')
121 FORMAT('# Hessian(i,i) / debye Angs^-2 amu^-1 ')
130 FORMAT('# Cubic(i,i,i) / hartree Angs^-3 amu^-3/2 ')
131 FORMAT('# Cubic(i,i,i) / debye Angs^-3 amu^-3/2 ')
140 FORMAT('# Quartic(i,i,i,i) / hartree Angs^-4 amu^-2 ')
141 FORMAT('# Quartic(i,i,i,i) / debye Angs^-4 amu^-2 ')

  !   >> >> 2 MR << <<

  tmp = ndof * (ndof - 1) / 2
  WRITE(8,201) title2MR

  WRITE(8,210)!# Hessian(i,j) / hartree Angs^-2 amu^-1
  DO i=2,ndof
     DO j=1,i-1
        WRITE(8,200)i,j, Hij(i,j)
     ENDDO
  END DO

  WRITE(8,220)!# Quartic(i,i,j,j) / hartree Angs^-4 amu^-2
  DO i=2,ndof
     DO j=1,i-1
        WRITE(8,200)i,j, Qiijj(i,j)
     END DO
  END DO

  WRITE(8,230)!# Cubic(i,i,j) / hartree Angs^-3 amu^-3/2
  DO i=2,ndof
     DO j=1,i-1
        WRITE(8,200)i,j, Ciij(i,j)
        WRITE(8,200)j,i, Ciij(j,i)
     END DO
  END DO

  WRITE(8,240)!# Quartic(i,i,i,j) / hartree Angs^-4 amu^-2
  DO i=2,ndof
     DO j=1,i-1
        WRITE(8,200)i,j, Qiiij(i,j)
        WRITE(8,200)j,i, Qiiij(j,i)
     END DO
  END DO

200 FORMAT(2i4,D20.10)
201 FORMAT('# 2MR ',a)
210 FORMAT('# Hessian(i,j) / hartree Angs^-2 amu^-1 ')
211 FORMAT('# Hessian(i,j) / debye Angs^-2 amu^-1 ')
220 FORMAT('# Quartic(i,i,j,j) / hartree Angs^-4 amu^-2 ')
221 FORMAT('# Quartic(i,i,j,j) / debye Angs^-4 amu^-2 ')
230 FORMAT('# Cubic(i,i,j) / hartree Angs^-3 amu^-3/2 ')
231 FORMAT('# Cubic(i,i,j) / debye Angs^-3 amu^-3/2 ')
240 FORMAT('# Quartic(i,i,i,j) / hartree Angs^-4 amu^-2 ')
241 FORMAT('# Quartic(i,i,i,j) / debye Angs^-4 amu^-2 ')
250 FORMAT(2i4,D20.10)

  !   >> >> 3 MR << <<
  tmp = ndof * (ndof - 1) * (ndof - 2) / 6
  WRITE(8,301) title3MR

  WRITE(8,310)!# Cubic(i,j,k) / hartree Angs^-3 amu^-3/2
  DO i=3,Ndof
     DO j=2,i-1
        DO k=1,j-1
           WRITE(8,300)i,j,k, Cijk(i,j,k)
        END DO
     END DO
  END DO

  WRITE(8,320)!# Quartic(i,i,j,k) / hartree Angs^-4 amu^-2
  DO i=3,Ndof
     DO j=2,i-1
        DO k=1,j-1
           WRITE(8,300)i,j,k, Qiijk(i,j,k)
           WRITE(8,300)j,k,i, Qiijk(j,k,i)
           WRITE(8,300)k,i,j, Qiijk(k,i,j)
        END DO
     END DO
  END DO

300 FORMAT(3i4,D20.10)
301 FORMAT('# 3MR ',a)
310 FORMAT('# Cubic(i,j,k) / hartree Angs^-3 amu^-3/2 ')
311 FORMAT('# Cubic(i,j,k) / debye Angs^-3 amu^-3/2 ')
320 FORMAT('# Quartic(i,i,j,k) / hartree Angs^-4 amu^-2 ')
321 FORMAT('# Quartic(i,i,j,k) / debye Angs^-4 amu^-2 ')

END SUBROUTINE write_QFF

SUBROUTINE make_quartic
  USE forceconstants
  USE global, ONLY: nMR
  IMPLICIT NONE
  Gi=0.0
  Ciii=0.0
  Qiiii=0.0
  IF(nMR .GT. 1) THEN
     Hij=0.0
     Ciij=0.0
     Qiiij=0.0
     IF(nMR .GT. 2) THEN
        Cijk=0.0
        Qiijk=0.0
     ENDIF
  ENDIF
END SUBROUTINE make_quartic

SUBROUTINE make_harmonic
  USE forceconstants
  USE global, ONLY: nMR
  IMPLICIT NONE
  Gi=0.0
  Ciii=0.0
  Qiiii=0.0
  IF(nMR .GT. 1) THEN
     Hij=0.0
     Ciij=0.0
     Qiiij=0.0
     Qiijj=0.0
     IF(nMR .GT. 2) THEN
        Cijk=0.0
        Qiijk=0.0
     ENDIF
  ENDIF
END SUBROUTINE make_harmonic

SUBROUTINE make_zero
  USE forceconstants
  USE global, ONLY: nMR
  IMPLICIT NONE
  Gi=0.0
  Hii=0.0
  Ciii=0.0
  Qiiii=0.0
  IF(nMR .GT. 1) THEN
     Hij=0.0
     Ciij=0.0
     Qiiij=0.0
     Qiijj=0.0
     IF(nMR .GT. 2) THEN
        Cijk=0.0
        Qiijk=0.0
     ENDIF
  ENDIF
END SUBROUTINE make_zero

SUBROUTINE put_coefs
  USE constants, ONLY: c_h2wn
  USE forceconstants
  USE global, ONLY: nMR
  IMPLICIT NONE
  ! Gi=Gi
  Hii=Hii * 0.5d0
  Ciii=Ciii  / 6.d0
  Qiiii=Qiiii  / 24.d0
  IF(nMR .GT. 1) THEN
     ! Hij=Hij
     Ciij=Ciij * 0.5d0
     Qiiij=Qiiij / 6.d0
     Qiijj=Qiijj  * 0.25d0
     IF(nMR .GT. 2) THEN
        !   Cijk=Cijk
        Qiijk=Qiijk  * 0.5d0
     ENDIF
  ENDIF
END SUBROUTINE put_coefs

SUBROUTINE convert2atomicunits
  USE constants
  USE forceconstants
  USE global, ONLY: nMR
  IMPLICIT NONE
  REAL*8:: converter
  converter=c_angs2bohr*dsqrt(c_amu2emass)
  Gi=Gi / converter
  Hii=Hii / converter/ converter
  Ciii=Ciii / converter/ converter/ converter
  Qiiii=Qiiii / converter/ converter/ converter/ converter
  IF(nMR .GT. 1) THEN
     Hij=Hij / converter/ converter
     Ciij=Ciij / converter/ converter/ converter
     Qiiij=Qiiij / converter/ converter/ converter/ converter
     Qiijj=Qiijj / converter/ converter/ converter/ converter
     IF(nMR .GT. 2) THEN
        Cijk=Cijk / converter/ converter/ converter
        Qiijk=Qiijk / converter/ converter/ converter/ converter
     ENDIF
  ENDIF
END SUBROUTINE convert2atomicunits

SUBROUTINE readinput
  USE global
  USE constants
  IMPLICIT NONE
  INTEGER,PARAMETER  :: n=50
  INTEGER::Nat,Nfree,m,MR,vmax_base,vmax(n)
  REAL*8 ::x(3*n),L(n*n),mass(n),omega(n)
  CHARACTER*2::Label(n)
  NAMELIST /mol/ x,L,mass,Label,Nat,omega,Nfree
  NAMELIST /vib/ MR, vmax, vmax_base, vscf, vci,vpt,vav,summary
  vmax=0
  MR=0
  vciorder=1
  vscf=.FALSE.
  vci=.FALSE.
  vpt=.FALSE.
  summary=.FALSE.
  REWIND(5)
  READ(5,mol)
  REWIND(5)
  READ(5,vib)
  nDOF=nfree
  nMR=MR
  IF(SUM(vmax) == 0 .AND. vmax_base /= 0)vmax=vmax_base

  ALLOCATE(maxbasis(Ndof),hrmfreq(Ndof),freq(Ndof))
  DO m=1,ndof
     hrmfreq(m)=omega(m)*c_wn2h
     maxbasis(m)=vmax(m)
     vciorder=vciorder*maxbasis(m)
  ENDDO
  PRINT*,hrmfreq*c_h2wn
  PRINT*,'nMR=',nMR
  PRINT*,'nDOF=',nDOF
  PRINT*,'hrmbasis=',maxbasis
  PRINT*,'vciorder=',vciorder
  PRINT*,'HO <0>=',SUM(hrmfreq)*0.5d0*c_h2wn
  freq=hrmfreq
END SUBROUTINE readinput

SUBROUTINE modeltest
  USE global
  USE forceconstants
  IMPLICIT NONE
  REAL*8::lambda,nu
  CALL make_zero()
  lambda=-0.1116d0
  nu=0.08414d0
  Hii(1)=0.5*(0.29375d0)
  Hii(2)=0.5*(2.12581d0)
  Ciij(2,1)=lambda
  Ciii(1)=lambda*nu
  hrmfreq(1)=dsqrt(0.29375d0)
  hrmfreq(2)=dsqrt(2.12581d0)
END SUBROUTINE modeltest

SUBROUTINE groundstate()
  USE global
  USE constants
  IMPLICIT NONE
  REAL*8 :: energy
  INTEGER::iter,state(nDOF)
  state=0
  CALL qVSCF(state,energy,iter)
  PRINT*,'Groud state energy is', energy*c_h2wn ,'after',iter,'iterations'
END SUBROUTINE groundstate

SUBROUTINE getVSCF(fundamentals)
  USE global
  USE constants
  IMPLICIT NONE
  REAL*8 :: energy,groundE
  REAL*8,OPTIONAL::fundamentals(0:nDOF)
  INTEGER::iter,state(nDOF),m
  state=0
  IF(VSCF)THEN
     PRINT*,"************************VSCF*****************"
     CALL myVSCF(state,energy,iter)
  ELSE
     PRINT*,"************************qVSCF*****************"
     CALL qVSCF(state,energy,iter)
  ENDIF
  IF (PRESENT(fundamentals))PRINT*, fundamentals
  groundE=energy
  PRINT*,'ground', 'is calculated after', iter, 'iterations for', state, groundE*c_h2wn
  IF (PRESENT(fundamentals))fundamentals(0)=energy
  PRINT*,'in cm-1'
  PRINT*,energy*c_h2wn
  DO m=1,nDOF
     state=0
     state(m)=1
     IF(VSCF)THEN
        CALL myVSCF(state,energy,iter)
     ELSE
        CALL qVSCF(state,energy,iter)
     ENDIF
     IF (PRESENT(fundamentals))fundamentals(m)=energy-groundE
     PRINT*,'fundamental', m, 'is calculated after', iter, 'iterations for', state, (energy-groundE)*c_h2wn
  ENDDO

  ! DO m=1,nDOF
  !    PRINT*,(fundamentals(m)-fundamentals(0))*c_h2wn
  ! ENDDO
END SUBROUTINE getVSCF

SUBROUTINE onemode(mode,energies,eigenfuns)
  USE global
  USE lapack
  USE yazar
  USE constants
  USE forceconstants
  USE integral
  IMPLICIT NONE
  REAL*8:: tmp,ham(maxbasis(mode),maxbasis(mode)),energies(MAXVAL(maxbasis))
  INTEGER:: mode,i,j,diff
  REAL*8, OPTIONAL,INTENT(out)::eigenfuns(maxbasis(mode),maxbasis(mode))
  Ham=0.d0
  PRINT*,'onemode'
  DO i=0, maxbasis(mode)-1
     tmp=0.d0
     tmp= tmp + (i + 0.5) * hrmfreq(mode) * 0.5;!kinetic term
     tmp= tmp + Hii(mode) * q_integral(i, mode, 2, 0);
     tmp= tmp + Qiiii(mode) * q_integral(i, mode, 4, 0);
     Ham(i+1,i+1)=tmp
     DO j=i+1,maxbasis(mode)-1!
        tmp=0.d0
        diff = j - i;
        IF (diff > 4) EXIT !break;
        IF (diff .EQ. 2) tmp= tmp-SQRT(float(j * (j - 1))) * hrmfreq(mode) * 0.25 !!kinetic term
        tmp= tmp + q_integral(j, mode, 1, diff) * gi(mode);
        tmp= tmp + q_integral(j, mode, 2, diff) * hii(mode);
        tmp= tmp + q_integral(j, mode, 3, diff) * ciii(mode);
        tmp= tmp + q_integral(j, mode, 4, diff) * qiiii(mode);
        Ham(i+1,j+1)=tmp
     ENDDO!!end for j
  ENDDO!!end for i

  CALL diag(maxbasis(mode),ham,energies)
  !if(present(eigenfuns))eigenfuns=ham !????????
  RETURN
END SUBROUTINE onemode

SUBROUTINE qonemr()
  USE global
  USE constants
  IMPLICIT NONE
  INTEGER :: mode
  REAL*8 :: groundE,energies(MAXVAL(maxbasis)),oneMRfunds(0:nDOF)
  REAL*8::HOfunds(0:nDOF),qVSCFfunds(0:nDOF)

  PRINT*,"********************qOne-MR"
  groundE=0.d0
  DO mode=1,nDOF
     CALL onemode(mode,energies)
     groundE=groundE+energies(1)
     oneMRfunds(mode)=energies(2)-energies(1)
     !print*,'mode ',mode,(oneMRfunds(mode))*c_h2wn
  ENDDO
  oneMRfunds(0)=groundE

  CALL getVSCF(qVSCFfunds)

  HOfunds(0)=SUM(hrmfreq)/2.d0
  HOfunds(1:ndof)=hrmfreq
  DO mode=0,nDOF
     PRINT*,'HO ',mode,(HOfunds(mode))*c_h2wn
     PRINT*,'1MR ',mode,(oneMRfunds(mode))*c_h2wn
     PRINT*,'qVSCF ',mode,(qVSCFfunds(mode))*c_h2wn
     PRINT*,'q1MR ',mode,(oneMRfunds(mode)-HOfunds(mode)+qVSCFfunds(mode))*c_h2wn
  ENDDO

END SUBROUTINE qonemr

SUBROUTINE myvci
  USE global
  USE forceconstants
  USE modvci
  USE integral
  USE lapack
  USE yazar
  USE constants
  IMPLICIT NONE
  INTEGER::i,j,m,m1,m2,m3,n,n1,n2,n3
  INTEGER::ketlabel(ndof),bralabel(ndof),diffvec(ndof),ndiffer,diff1,diff2,diff3,diff4
  REAL*8:: tmp,ham(vciorder,vciorder),energies(vciorder)
  ham=0.d0
  DO i=0,vciorder-1
     tmp=0.d0
     CALL statelabel(i,ketlabel)
     DO m=1,ndof
        n = ketLabel(m);
        tmp=tmp+ hii(m) * q_integral(n, m, 2, 0) + qiiii(m) * q_integral(n, m, 4, 0);
        tmp=tmp+ (ketLabel(m) + 0.5) * hrmfreq(m) * 0.5;
     ENDDO
     IF (nMR > 1)THEN
        DO m1=1,nDOF
           n1 = ketLabel(m1);
           DO m2=m1+1, nDOF
              n2 = ketLabel(m2);
              tmp=tmp+ qiijj(m1,m2) * q_integral(n1, m1, 2, 0) * q_integral(n2, m2, 2, 0);
           ENDDO!m2
        ENDDO !m1
     ENDIF !nMR>1
     Ham(i+1,i+1)=tmp
     DO j=i+1,vciorder-1
        tmp=0.d0
        CALL stateLabel(j,bralabel);
        CALL fundiffer(bralabel,ketlabel,diffvec,ndiffer)
        IF (nDiffer > 3) THEN
           CYCLE
        ELSE IF (nDiffer == 1)THEN
           diff1 = abs(braLabel(diffVec(1)) - ketLabel(diffVec(1)));
           n1 = MAX(braLabel(diffVec(1)), ketLabel(diffVec(1)));
           tmp=tmp+ q_integral(n1, diffVec(1), 1, diff1) * gi(diffVec(1));
           tmp=tmp+ q_integral(n1, diffVec(1), 2, diff1) * hii(diffVec(1));
           tmp=tmp+ q_integral(n1, diffVec(1), 3, diff1) * ciii(diffVec(1));
           tmp=tmp+ q_integral(n1, diffVec(1), 4, diff1) * qiiii(diffVec(1));
           IF (diff1 == 2) tmp=tmp - DSQRT(Dfloat(n1 * (n1 - 1))) * hrmfreq(diffVec(1)) * 0.25!kinetic term
           IF (nMR > 1)THEN
              DO m=1,ndof
                 IF (m == diffVec(1))CYCLE
                 n2 = MAX(braLabel(m), ketLabel(m));
                 diff2 = abs(braLabel(m) - ketLabel(m));
                 tmp=tmp+ qiijj(m,diffVec(1)) * q_integral(n2, m, 2, diff2) * q_integral(n1,diffVec(1), 2, diff1);
                 tmp=tmp+ ciij(m,diffVec(1)) * q_integral(n2, m, 2, diff2) * q_integral(n1,diffVec(1), 1, diff1);
              ENDDO !m1
           ENDIF!nmr>1
        ELSE IF ((nDiffer == 2) .AND. (nMR > 1)) THEN
           m1 = diffVec(1);
           m2 = diffVec(2);
           diff1 = abs(braLabel(m1) - ketLabel(m1));
           diff2 = abs(braLabel(m2) - ketLabel(m2));
           n1 = MAX(braLabel(m1), ketLabel(m1));
           n2 = MAX(braLabel(m2), ketLabel(m2));
           tmp=tmp+ hij(m1,m2) * q_integral(n1, m1, 1, diff1) * q_integral(n2, m2, 1, diff2);
           tmp=tmp+ qiijj(m1,m2) * q_integral(n1, m1, 2, diff1)* q_integral(n2, m2, 2, diff2);
           tmp=tmp+ ciij(m1,m2) * q_integral(n1, m1, 2, diff1) * q_integral(n2, m2, 1, diff2);
           tmp=tmp+ qiiij(m1,m2) * q_integral(n1, m1, 3, diff1)* q_integral(n2, m2, 1, diff2);
           tmp=tmp+ ciij(m2,m1) * q_integral(n1, m1, 1, diff1) * q_integral(n2, m2, 2, diff2);
           tmp=tmp+ qiiij(m2,m1) * q_integral(n1, m1, 1, diff1)* q_integral(n2, m2, 3, diff2);
           IF (nMR > 2)THEN
              DO m=1,ndof
                 IF (m == m1 .OR. m == m2)CYCLE;
                 n = braLabel(m);
                 tmp=tmp+ qiijk(m,m1,m2) * q_integral(n, m, 2, 0) * q_integral(n1, m1, 1, diff1)&
                      * q_integral(n2, m2, 1, diff2);
              ENDDO
           ENDIF

        ELSE IF ((nDiffer == 3) .AND. (nMR > 2))THEN
           m1 = diffVec(1);
           m2 = diffVec(2);
           m3 = diffVec(3);
           diff1 = abs(braLabel(m1) - ketLabel(m1));
           diff2 = abs(braLabel(m2) - ketLabel(m2));
           diff3 = abs(braLabel(m3) - ketLabel(m3));
           n1 = MAX(braLabel(m1), ketLabel(m1));
           n2 = MAX(braLabel(m2), ketLabel(m2));
           n3 = MAX(braLabel(m3), ketLabel(m3));
           tmp=tmp+ cijk(m1,m2,m3) * q_integral(n1, m1, 1, diff1)* q_integral(n2, m2, 1, diff2)&
                * q_integral(n3, m3, 1, diff3);
           tmp=tmp+ qiijk(m1,m2,m3) * q_integral(n1, m1, 2, diff1) * q_integral(n2, m2, 1,diff2)&
                * q_integral(n3, m3, 1, diff3);
           tmp=tmp+ qiijk(m2,m1,m3) * q_integral(n2, m2, 2, diff2) * q_integral(n1, m1, 1,diff1)&
                * q_integral(n3, m3, 1, diff3);
           tmp=tmp+ qiijk(m3,m2,m1) * q_integral(n3, m3, 2, diff3) * q_integral(n2, m2, 1,diff2)&
                * q_integral(n1, m1, 1, diff1);
        ENDIF
        ham(i+1,j+1)=tmp
     ENDDO!end for j
  ENDDO!end for i
!  CALL writematrix(vciorder,ham*c_h2wn)
  CALL diag(vciorder,ham,energies)
  PRINT*,'vci energies'
  PRint*,energies(1)*c_h2wn
  do m=1,ndof*3
    print*, (energies(m+1)-energies(1))*c_h2wn
  enddo
END SUBROUTINE myvci

PROGRAM main
  USE global
  USE qVSCFintegrals
  USE forceconstants
  USE constants
  IMPLICIT NONE

  CALL readinput()
  CALL init_qff()
  CALL read_qff()
  ! CALL make_quartic()
  ! CALL make_harmonic()
  !CALL make_zero
  CALL convert2atomicunits()
  ! CALL write_qff()
  ! CALL groundstate()
  CALL put_coefs()!only for VSCF
  CALL qonemr()
  ! CALL getVSCF()
  CALL myvci()
  STOP
END PROGRAM main
